clear
load ../data/FLAC_results/ld_profiles.mat
peakfigs=1; %Display figures of peak locations or not
F.Z=squeeze(Z(:,1,:)); %bathymetry
F.X=X(1,:); %distance from ridge axis
F.U=linspace(3.0,30.3,40)*1e-10; %spreading half-rates, m/s

%Convert to time & trim to 780 ky
F.t1=0:780; F.t2=-780:0; F.t=-780:780;
for i=1:length(F.Z(:,1))
    if i>40; j=i-40; else; j=i; end;
    F.Zt1(i,:)=interp1(3.171e-8*F.X/F.U(j),F.Z(i,:),F.t1);
    F.Zt2(i,:)=interp1(3.171e-8*F.X/F.U(j),F.Z(i,:),F.t2);
end 

%run to display
n=47;

%Calculate spacing
p=0.025; %minimum peak prominence

%Right side
for i=1:length(F.Zt1(:,1))
    Z=F.Zt1(i,:);
    [P,k]=findpeaks(Z,F.t1,'MinPeakProminence',p);
    
    %add missed faults manually for slow spreading ridges
    if i==5; k(2)=774.4; P(2)=nan; end;
    if i==7; k(3)=747.1; P(3)=nan; end;
    if i==45; k(3)=774.4; P(3)=nan; end;
    
    %fault spacing in time
    a=diff(k);
    
    %include faults closest to axis for slow spreading ridges
    if (i<8) | (i>40 & i<48)
        F.pdist1(i)={a(1:end)};
    else
        F.pdist1(i)={a(2:end)};
    end
    
    if length(k)>0; F.in1(i)=k(1); end;
    
    if peakfigs==1 && i==n % i = profile number you want to display
    figure('Position',[500,500,1500,400]);
    set(gca,'FontSize',16);
    hold on; 
    plot(k,P,'*','MarkerSize',10);
    plot(X(1,:)./(F.U(rem(i,40))*3.154e7),F.Z(i,:),'LineWidth',2,'Color',[0.8 0 0]);
    xlim([0,780]);
    xlabel('Time (kyr)');
    ylabel('Height (km)');
    hold off;
    end 
    
end 

%Left side
for i=1:length(F.Zt2(:,1))
    Z=F.Zt2(i,:);
    [P,k]=findpeaks(Z,F.t2,'MinPeakProminence',p);
   
    %add missed faults manually for slow spreading ridges
    %none to report yet

    a=diff(k);

    %include fault closest to axis for slowest spreading ridges
    if (i<8) | (i>40 & i<48)
        F.pdist2(i)={a(1:end)};  
    else
        F.pdist2(i)={a(1:end-1)}; %Note removal of closest to axis 
    end
    
    if length(k)>0; F.in2(i)=k(end); end;
    
    if peakfigs==1 && i==n % i = profile number you want to display
    figure('Position',[500,500,1500,400]);
    set(gca,'FontSize',16);
    hold on; 
    plot(k,P,'*','MarkerSize',10);
    plot(X(1,:)./(F.U(rem(i,40))*3.154e7),F.Z(i,:),'LineWidth',2,'Color',[0.8 0 0]);
    xlim([-780,0]);
    xlabel('Time (kyr)');
    ylabel('Height (km)');
    hold off;
    end 
end 

%Combine
for i=1:length(F.Z(:,1))
    combine=[cell2mat(F.pdist2(i)), cell2mat(F.pdist1(i))];
    if sum(isnan(combine))==0
    F.pdist(i)=mean(combine);
    F.prange(i)=std(combine);
    else
    plf=find(~isnan(combine));
    F.pdist(i)=combine(plf);
    F.prange(i)=nan;
    end
end 

%Split into forcing periods
F.pr41=F.pdist(1:40);
F.pr100=F.pdist(41:80);

F.U1=F.U*3.154e7;
save('spacingld.mat','F');

